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ABSTRACT 


Fish of the superfamily Cobitoidea sensu stricto 
(namely loaches) exhibit extremely high diversity of 
color patterns, but so far little is known about their 
evolutionary mechanism. Melanocortin 1 receptor 
gene (MC1R) plays an important role during the 
synthesis of melanin and formation of animal body 
color patterns. In this study, we amplified and 
sequenced the partial MC1R gene for 44 loach 
individuals representing 31 species of four families. 
Phylogenetic analyses yielded a topology congruent 
with previous studies using multiple nuclear loci, 
showing that each of the four families was 
monophyletic with sister relationships of Botiidae+ 
(Cobitidae+(Balitoridae+Nemacheilidae)). Gene 
evolutionary analyses indicated that MC1R in 
loaches was under purifying selection pressure, with 
various sites having different dy/ds values. Both 
Botiidae and Cobitidae had lower dy/ds values than 
those of background lineages, suggesting their 
evolution might be strongly affected by purifying 
selection pressure. For Balitoridae and 
Nemacheilidae, both had larger dyds values than 
those of background lineages, suggesting they had a 
faster evolutionary rate under more relaxed selection 
pressure. Consequently we inferred that the 
relatively stable color patterns in Botiidae and 
Cobitidae might result from the strong purifying 
selection pressure on the MC1R gene, whereas the 
complicated and diverse color patterns in Balitoridae 
and Nemacheilidae might be associated with the 
relaxed selection pressure. Given the easy 
experimental procedure for the partial MC1R gene 
and its excellent performance in reconstructing 
phylogeny, we suggest this gene could be used as a 
good molecular marker for the phylogenetic study of 
fish species. 
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INTRODUCTION 


The superfamily Cobitoidea is a group of small- to medium- 
Sized benthic fish, composed of approximately 2896 of species 
of the order Cypriniformes, which is the largest group of 
freshwater fish in the world (Nelson et al., 2016). Depending on 
different authors, Cobitoidea includes variable families. Bohlen 
& Slechtová (2009) and Chen et al. (2009) congruently 
recognized the genus Ellopostoma as a distinct new family 
Ellopostomatidae, and proposed that Cobitoidea is composed 
of eight families (Catostomidae, Gyrinocheilidae, Botiidae, 
Vaillantellidae, Cobitidae, Ellopostomatidae, Nemacheilidae and 
Balitoridae). Kottelat (2012) raised genera Serpenticobitis and 
Barbucca to family rank, and established Serpenticobitidae and 
Barbuccidae. Both Serpenticobitis and Barbucca have been 
formerly included in Balitoridae (Bohlen & Slechtová, 2009; 
Slechtová et al., 2007), and consist of three and two species, 
respectively. While the most recent version of Fishes of the 
World (Nelson et al., 2016) follows Kottelat (2012), without 
more convincing evidence, we herein adopt the classification of 
Cobitoidea suggested by Slechtová et al. (2007), Bohlen & 
Slechtová (2009), and Chen et al. (2009). 

Cobitoidea sensu stricto, which excludes Gyrinocheilidae 
(algae eaters) and Catostomidae (suckers), is commonly known 
as the loaches (Nelson et al., 2016). This group of fish exhibit 
extremely high color pattern diversity, resulting in the observed 
abundant biodiversity of loaches. This makes them ideal 
subjects for studying adaptive evolution to environments, 
camouflage, kin recognition, and mate choice. Their beautiful 
color patterns also attract numerous fish enthusiasts, making 
them one of the most important groups of aquarium fish. A 
number of community websites about loaches have been 
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established by fish hobbyists, with the most detailed site 
(Loaches Online, http: //www.loaches.com/) containing 
comprehensive knowledge on their taxonomy, form, disease, 
and care. Additionally, some loach species are considered as 
important economic fish, including Misgurnus anguillicaudatus 
(mud loach), which is widely farmed in China (Zhang et al., 
2012). Therefore, loaches play significant roles in both 
evolutionary biology and social economy. To date, however, the 
phylogenetic relationships of loaches still remain controversial 
(Liu et al., 2012; Mayden et al., 2009; Slechtova et al., 2007; 
Tang et al., 2006). As for their color patterns, almost no 
investigations have been performed, and little is known about 
their evolutionary mechanisms. 

Among the many genes related to color patterns, 
melanocortin 1 receptor gene (MC1R) is one of the most widely 
studied, belonging to a family of G-protein-coupled receptors 
involved in various physiological processes in vertebrates and 
playing an important role in the synthesis of melanin and 
formation of animal body color patterns (Selz et al., 2007). The 
MC1R gene has been extensively investigated in numerous 
animal species, including reptiles (Corso et al., 2012; Cox et al., 
2013; Herczeg et al., 2010; Nunes et al., 2011; Rosenblum et 
al, 2004), birds (Baiao & Parker, 2012; Mundy, 2005), 
mammals (Ayoub et al., 2009; Lu & Zhang, 2001; Pérez et al., 
2013), and fish (Bar et al., 2013; Gross et al., 2009; Henning et 
al., 2010; Selz et al., 2007). Almost all fish species studied thus 
far have shown that the MC1R gene is conserved as a single 
copy and single exon protein-coding gene with a 966 bp long 
open reading frame (ORF), including seven transmembrane 
domains (Bar et al., 2013; Braasch et al., 2008; Henning et al., 
2010; Hofreiter & Schóneberg, 2010; Nunes et al., 2011), 
suggesting that this gene should be a good molecular marker 


Table 1 Samples used in the present study 


for reconstructing fish phylogeny. So far, however, only a few 
studies have applied this gene to mammal phylogenies (Ayoub 
et al., 2009; Pérez et al., 2013), and no report exists for fish 
phylogenies. The color patterns of loaches are highly diverse, 
which make them good subjects to study the evolution of the 
MC1R gene. In the present study, we amplified and sequenced 
the partial sequence of the MC7R gene in loaches, and 
analyzed their characteristics. Using this gene as a molecular 
marker, the phylogeny of loaches was reconstructed, and the 
evolution of the MC1R gene was then analyzed. 


MATERIALS AND METHODS 


Sample collection 

Because species of the families Ellopostomatidae and 
Vaillantellidae are distributed in Borneo, Sumatra, and Thailand, 
and only consist of two and three species, respectively, 
samples can be difficult to obtain. In this study, four families of 
loaches were involved. In total, 52 individuals representing 34 
species in Cypriniformes were selected for analyses of MC1R 
gene evolution, which included 44 individuals belonging to 31 
loach species and several outgroups, including one 
Gyrinocheilus aymonieri, one Hypophthalmichthys molitrix, and 
six Danio rerio individuals. To guarantee the accuracy of the 
targeted sequences, the cDNA sequences of four D. rerio 
individuals were downloaded from GenBank. In addition, we 
sequenced two D. rerio individuals, and all 46 other sequences 
in the present study. All samples were collected from the main 
drainages in China and aquarium markets, preserved in 95% 
ethanol, and deposited in the Institute of Hydrobiology, Chinese 
Academy of Sciences. Detailed information of all samples is 
listed in Table 1, with GenBank accession numbers included. 





Family Species and individuals Locality GenBank accession No. 
Botiidae Sinibotia reevesae Chishui, Guizhou Province KX120160 
Yasuhikotakia modesta Market in Wuhan KX120161 
Leptobotia elongata Market in Wuhan KX120164 
Leptobotia pellegrini Quzhou, Zhejiang Province KX120167 
Leptobotia microphthalma 1 Tucheng, Guizhou Province KX120162 
Leptobotia microphthalma 2 Tucheng, Guizhou Province KX120163 
Leptobotia taeniops 1 Yuanjiang, Hunan Province KX120165 
Leptobotia taeniops 2 Yuanjiang, Hunan Province KX120166 
Parabotia maculosus 1? Jianou, Fujian Province KX120168 
Parabotia maculosus 2° Jianou, Fujian Province KX120169 
Parabotia fasciatus Tucheng, Guizhou Province KX120170 
Cobitidae Cobitis sinensis 1 Wuyishan, Fujian Province KX120157 
Cobitis sinensis 2 Wuyishan, Fujian Province KX120158 
Pangio kuhlii Market in Wuhan KX120159 
Nemacheilidae Homatula wujiangensis Chishui, Guizhou Province KX120149 
Homatula variegata 4n Chishui, Guizhou Province KX120146 
Homatula variegata 2 Baoji, Shaanxi Province KX120148 
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Family Species and individuals 


Continued 





Nemacheilidae Homatula variegata 3° 
Schistura fasciolata 
Schistura incerta 
Traccatichthys taeniatus 1 
Traccatichthys taeniatus 2 
Triplophysa bleekeri 
Triplophysa sp. 1 
Triplophysa sp. 2 
Balitoridae Erromyzon kalotaenia 1° 
Erromyzon kalotaenia 2° 
Vanmanenia stenosoma 
Vanmanenia pingchowensis 
Vanmanenia lineata 
Vanmanenia hainanensis 
Vanmanenia caldwelli 1° 
Vanmanenia caldwelli 2° 
Formosania stigmata 
Formosania chenyiyui 1 
Formosania chenyiyui 2 
Pseudogastromyzon fasciatus 1 
Pseudogastromyzon fasciatus 2 
Pseudogastromyzon fangi 1° 
Pseudogastromyzon fangi 2° 
Beaufortia kweichowensis 
Jinshaia abbreviata ' 

Jinshaia sinensis ' 
Sinogastromyzon szechuanensis 
Outgroups Gyrinocheilus aymonieri 
Hypophthalmichthys molitrix 
Danio rerio 1 

Danio rerio 2 

Danio rerio 3 

Danio rerio 4 

Danio rerio 5 


Danio rerio 6 


Individuals with the same superscript letters share the same haplotype. 


DNA extraction, PCR amplification, and sequencing 

Total DNA was extracted from muscles following the salt- 
extraction procedure of Aljanabi & Martinez (1997), with some 
modification as per Tang et al. (2008). Five pairs of primers 
were newly designed according to the conserved regions of the 
aligned teleost MC1R gene sequences, and the lengths of the 


Locality GenBank Accession No. 
Baoji, Shaanxi Province KX120147 
Xichang, Sichuan Prov. KX120150 
Changting, Fujian Province KX120151 
Zhongshan, Guangzhou Province KX120156 
Zhongshan, Guangzhou Province KX120155 
Chishui, Guizhou Province KX120154 
Baoji, Shaanxi Province KX120152 
Baoji, Shaanxi Province KX120153 
Laibin, Guangxi Province KX120133 
Laibin, Guangxi Province KX120134 
Quzhou, Zhejiang Province KX120132 
Yuanjiang, Hunan Province KX120127 
Laibin, Guangxi Province KX120129 
Yuedong, Hainan Province KX120128 
Shaowu, Fujian Province KX120138 
Shaowu, Fujian Province KX120139 
Huangken, Fujian Province KX120142 
Changting, Fujian Province KX120130 
Changting, Fujian Province KX120131 
Guangze, Fujian Province KX120135 
Longyan, Fujian Province KX120136 
Laibin, Guangxi Province KX120140 
Longyan, Fujian Province KX120141 
Market in Wuhan KX120137 
Panzhihua, Sichuan Province KX120144 
Panzhihua, Sichuan Province KX120143 
Chishui, Guizhou Province KX120145 
Market in Wuhan KX120173 
Market in Wuhan KX120174 
Market in Wuhan KX120172 
Market in Wuhan KX120171 
AY161847* 
NM180970* 
BC162848* 
BC162836* 


*: sequences downloaded from GenBank. 


target segments ranged from 650 bp to 850 bp. However, only 
by using the primer set for the shortest target segment could we 
successfully obtain sequences for all samples. Thus, the primer 
set used was inMC1RF (5-AGCGTCAGYAAYGTGGTGGAGA- 
3) and inMC1RH (5-CGGTTCTGTACCTGCACAT-3?. 

The polymerase chain reaction (PCR) mixtures contained 3 uL of 
10x taq Buffer, 0.9 uL of dNTPs (10 mmol/L), 0.75 uL of each 
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primer (10 mmol/L), 0.15 uL of Taq polymerase, 1 uL of DNA 
template, and sterile distilled H2O to a final volume of 30 pL. 
The PCR thermal cycle profile was as follows: initial denaturation 
step at 94 °C for 5 min, followed by 35 cycles at 94 °C for 30 sec, 
56 °C for 45 sec, 72 °C for 1 min, and a final extension at 72 °C 
for 8 min. The amplified fragments were purified and sequenced 
by Sangon Biotech (Shanghai) Co., Ltd. To confirm the 
accuracy of sequencing, sequencing reactions were performed 
from both ends of each fragment using the same primers as 
PCR amplification. 


Sequence variation and phylogenetic analyses 

Multiple alignments of sequences were performed using 
ClustalX 2.1 (Larkin et al., 2007; Thompson et al., 1997) 
alignment editor with the default settings, with all aligned 
sequences then translated into amino acid residues in MEGA 
6.0 (Tamura et al., 2013) to test for sequencing mistakes. The 
alignments were verified by eye in SEAVIEW (Galtier et al., 
1996). Saturation of genes and codon positions were assessed 
using the ‘transitions and transversions vs. divergence’ graphic 
function in DAMBE 5 (Xia, 2013). Nucleotide compositions and 
pairwise genetic distances were calculated based on the 
Kimura 2-parameter model with standard errors estimated 
using 1 000 bootstrap replicates. Molecular phylogenetic 
relationships were estimated using Bayesian inference (BI) and 
maximum likelihood (ML) methods. The optimal model of 
nucleotide evolution for BI analyses was identified using the 
Bayesian Information Criterion (BIC) as estimated in jModelTest 
2.1.7 (Darriba et al., 2012; Guindon & Gascuel, 2002). The 
TPM1uf+G model was selected as the best model, so we set 
the parameter nst=6 with rates=gamma when running the BI 
analyses. 

Bayesian inference (Bl) was carried out using MrBayes 3.1.2 
(Ronquist & Huelsenbeck, 2003). Two independent analyses 
with four simultaneous Markov chains were run for 4 000 000 
generations, sampling every 1 000 generations, with a total of 4 
001 trees each. We discarded all samples obtained during the 
first 1 000 000 generations as "burn-in", and a 5096 majority- 
rule consensus tree with posterior probability values for each 
node was obtained from the remaining 3 001 trees (Figure 1). 

Maximum likelihood (ML) analysis was conducted using 
RAxML 8.1.13 (Stamatakis, 2014). We searched for the best 
scoring ML tree using a general time-reversible nucleotide 
model (GTR) with gamma-distributed rate variation among sites 
(G) and invariable sites (Il). A rapid bootstrap analysis with 1 
000 replicates was used to assess the relative robustness of 
node support. 


Analyses of MC1R gene evolution 

The w ratio (defined as the ratio of non-synonymous vs. 
synonymous substitution rates, i.e., w=dy/ds) is a measure of 
natural selection acting on a protein. Simplistically, values of 
w<1, 21, and >1 means negative purifying selection, neutral 
evolution, and positive selection, respectively (Yang, 2007). Site 
models and branch models were separately adopted to 
estimate whether the MC1R gene in the present study had 
experienced different selection pressures along different sites 
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and different lineages. The CODEML program within the PAML 
package was used to assess parameters in models of 
sequence evolution and to test relevant hypotheses (Yang, 
2007). In the following analyses, the outgroups were 
excluded, and the input tree was an unrooted tree with 
trifurcation at the root. 

For site models, we examined three pairwise codon-based 
substitution models to assess w values for all codon sites: MO 
(one-ratio) vs. M3 (discrete w), M1a (nearly neutral) vs. M2a 
(positive selection), and M7 (8 distribution) vs. M8 (B distribution 
and a fraction of sites with w>1). Likelihood ratio tests (LRTs) 
were performed to compare the fit of two pairwise models. It is 
assumed that twice the log likelihood difference between 
nested models (2AInL) follows a X^ distribution with the number 
of degrees of freedom equal to the difference in the number of 
free parameters (Whelan & Goldman, 1999). When LRTs 
indicate positive selection, the Bayesian empirical Bayesian 
(BEB) method (Yang et al., 2005) can be used to calculate the 
posterior probabilities of positively selected codons. 

For branch models, based on the phylogenetic tree yielded in 
the present study (Figure 1), we divided loaches into four 
lineages corresponding to four families (Balitoridae, 
Nemacheilidae, Cobitidae, and Botiidae). The following steps 
were performed: (1) w values of the MC1R gene dataset were 
estimated for all lineages using a one-ratio model (all lineages 
share a common w value); (2) each lineage was separately 
treated as a foreground lineage compared with the remaining 
lineages (background lineages), and  two-ratio models 
(foreground lineage and background lineages have different w 
values) were independently conducted for each lineage; (3) four 
pairs of LRTs were examined between the one-ratio and two- 
ratio models for four lineages. 


RESULTS 


Characteristics of MC1R gene sequences 

After alignment, 591 bp of the MC1R gene sequences was 
used for analysis. No indels or deletions were detected. In total, 
43 haplotypes were identified from 52 individuals. For the 44 
loach individuals, 38 haplotypes were recognized. Except for 
Jinshaia abbreviata and J. sinensis sharing the same haplotype, 
all other shared haplotypes occurred among different 
individuals of the same species. The average nucleotide 
composition for all species was A=19.5%, T=28.5%, G=22.1%, 
and C=29.9%. The content of A+T (48.0%) was lower than that 
of C+G (52.0%). Strong compositional biases against A existed 
at the third position (only 7.6%). For the 591 bp sequences, 226 
sites were variable, of which 194 were parsimony informative. 
The average ratio of transitions/transversions (Ti/Tv) was 2.084. 
Plots of the absolute transitions and transversions versus TrN 
distance showed that the transitions and transversions for all 
three codon positions did not reach saturation (not shown). 


Phylogenetic relationships of Cobitoidea 

Maximum likelihood and BI analyses yielded almost the same 
topologies, except that the lineage (representing one botiid 
subfamily, Botiinae) formed by Sinibotia reevesae and 
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Figure 1 ML tree estimated using RaxML for all individuals based on partial MC1R gene sequences 
Values at the nodes correspond to support values for ML/BI methods (only values >50 or 0.50 shown). *: a support value of 100 or 1.00; -: values of «50 
or 0.50. Species names are followed by individual codes as in Table 1. 
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Yasuhikotakia modesta had a different position. In the ML tree, 
this lineage clustered with botiid subfamily Leptobotiinae, 
formed by Leptobotia and Parabotia species, and together they 
formed the family Botiidae. Whereas in the BI tree, Botiinae was 
located in the basal position of all loaches, but with a very low 
posterior probability (only 0.44). The ML tree with node support 
values of both methods labeled is shown in Figure 1. 

The phylogenetic tree indicated that each of the four loach 
families was monophyletic and had high support values 
(bootstrap value (BP) higher than 94 in ML, and all posterior 
probability (PP)=1.00 in Bl), except for Botiidae. The sister 
relationships among the four families was Botiidae* (Cobitidae- 
(Balitoridae+Nemacheilidae)). Balitoridae and Nemacheilidae 


formed a sister group with medium support values (BP=91 in ML, 
and PP=0.89 in BI), located at the crown of the phylogenetic 
tree. The family Cobitidae was at an intermediate position. 
Within Botiidae, Botiinae was distantly related to the 
Leptobotiinae with a genetic distance of 0.092+0.013. For 
balitorid species, subfamilies Balitorinae and Gastromyzontinae 
were both monophyletic and formed monophyletic Balitoridae. 
Within Gastromyzontinae, however, the species of genera 
Vanmanenia, Formosania, and Pseudogastromyzon were not 
monophyletic. The genetic distance between these two 
subfamilies was not very high (0.073+40.010). The pairwise 
mean genetic distances of the MC1R gene among the loach 
groups are listed in Table 2. 


Table 2 Pairwise mean genetic distance (lower triangle) and standard error (upper triangle) based on the Kimura 2-parameter model for 


loach groups and outgroups 








Group code 
Family or subfamily Group code 

1 2 3 4 5 6 7 
Gastromyzontinae 1 0.010 0.012 0.012 0.014 0.013 0.018 
Balitorinae 2 0.073 0.011 0.012 0.014 0.012 0.017 
Nemacheilidae 3 0.114 0.093 0.013 0.015 0.012 0.017 
Cobitidae 4 0.114 0.090 0.121 0.013 0.012 0.016 
Botiinae 5 0.123 0.120 0.142 0.118 0.011 0.016 
Leptobotiinae 6 0.116 0.093 0.104 0.101 0.092 0.015 
Outgroups 7 0.190 0.174 0.186 0.175 0.167 0.150 


Bold indicates the genetic distance between two subfamilies of the same family. 


Evolution of the MC1R gene in Cobitoidea 

Using the topology yielded by the ML method, we analyzed the 
evolution of MC1R in Cobitoidea. Analysis based on the one- 
ratio model MO showed that the average w value of the MC1R 
gene in all loaches was 0.034 9, much less than 1, suggesting 
that the average selection pressure on this gene was negative 
or purifying. Using LRT analyses, comparisons of the pairwise 
models indicated that the M3 model (having discrete w values 
among sites; -InL=3 295.024 728) was significantly better than 
the MO model (having the same w value of all sites; -/nL=3 
366.540 405) (P=0.000) at estimating w values among sites; 


however, there was no obvious difference between M1a 
(nearly neutral; -/nL-3 301.622 991) and M2a (positive 
selection; -/nL-3 301.622 991) models (P=1.000), suggesting 
that no amino acid site was obviously under positive 
selection, although different sites had variable 
evolutionary rates. However, comparison also indicated 
that the M8 model (-/nL-3 297.026 100) was significantly 
better than the M7 model (-/nL-3 304.204 399) (P=0.001), 
and the BEB procedure identified amino acid sites 12Q, 20Q, 
and 176H as possibly being under positive selection, though 
each posterior probability (Pr) was less than 95% (59.2%, 
73.096, and 87.396, respectively) (Table 3). 


Table 3 Likelihood ratio tests (LRTs) between selected CODEML codon substitution models for the MC1R gene of Cobitoidea 


2AInL 
143.033 2 


Compared models InL 
M3-MO M3: -3 295.024 728 


MO: -3 366.540 405 


M2a-M1a M2a: -3 301.622 991 0 
Mia: -3 301.622 991 
M8-M7 M8: -3 297.026 100 14.356 6 


MT: -3 304.204 399 


For variation of the w values among the loach lineages, 
pairwise comparisons of LRTs between so called "foreground" 
and "background" lineages showed that two pairs were 
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df P value Positive selected sites (posterior probability) 
4 0.000 

2 1.000 

2 0.001 12Q (59.2%), 20Q (73.0%), 


176H (87.3%) 


detected with significant variations. The mean w value of 
Balitoridae (0.055 1) was significantly higher than that of the 
other three families (denoted as non-balitorids, 0.026 0) 


(P=0.006), whereas the mean w value of Cobitidae (0.002 8) 
was significantly lower than that of the non-cobitids (0.040 5) 
(P=0.00). For Botiidae, its mean w value was lower (0.025 2) 
than that of the non-botiids (0.039 8) (P=0.096). The mean w 


value of the nemacheilids (0.046 0) was slightly higher than 
that of the non-nemacheilids (0.032 1) (P=0.211). The 
pairwise comparison results for the four lineages are listed in 
Table 4. 


Table 4 LRTs performed to detect heterogeneous selection regimes among lineages for the MC1R gene 








Branch-specific model Parameter estimates df InL 2AInL P value 

All lineages 

MO (one-ratio) w=0.034 9 87  -3 366.540 405 

Balitoridae* 

two-ratio vs. one-ratio Wo=0.026 0 w1=0.055 1 88 | -3 362.719 906 7.640 998 0.006<0.01** 


Nemacheilidae? 
two-ratio vs. one-ratio Wo=0.032 1 w1=0.046 0 
Cobitidae? 

two-ratio vs. one-ratio Wo=0.040 5 w1=0.002 8 
Botiidae? 


two-ratio vs. one-ratio (970.039 8 (470.025 2 


88 | -3365.759 399 1.562 012 0.211 


88 . —-3358.153 163 16.774 48 0.00«0.01** 


88 | -3365.156 736 2.767 338 0.096 


a: foreground lineage; w4: value of dw/ds for foreground lineages; wo: value of dw/ds for background lineages. 


DISCUSSION 


Evolutionary variation of MC1R among sites and loach 
lineages 

Presently, dozens of studies have indicated that the MC1R 
gene is highly conserved among vertebrates, and variations are 
closely associated with pigmentation differences (Hubbard et al., 
2010; Mundy, 2005; Nunes et al., 2011; Rosenblum et al., 2004), 
with only a few exceptions (Cheviron et al., 2006; Dorn et al., 
2011). Such findings are important in clarifying the origin of new 
species or local adaptation within species (Hubbard et al., 
2010). Compared with other vertebrates, teleost fish have more 
diverse color patterns (Parichy, 2003) Among numerous 
cypriniform fish, loaches are one of the most important groups 
with colorful pigment patterns. 

The present results showed that in loaches the MC1R gene 
was under negative selection pressure overall (average value 
of w=0.034 9«1). Analyses based on site models indicated that 
the M3 model was significantly better than the MO model 
(P=0.000), suggesting that various sites of the MC1R gene had 
disparate selection pressure. However, no amino acid site was 
obviously under positive selection since no significant difference 
existed between the Mia (nearly neutral) and M2a (positive 
selection) models (P=1.000). In contrast to the M1a-M2a 
comparison, the results of the M7-M8 comparison were 
significant and model M8 detected three positive amino acid 
sites, though with low posterior probabilities. Yang (2007) 
determined that the w estimate under M8 might often produce 
false positives and the M1a-M2a comparison could be more 
robust. Therefore, given the results of the M1a-M2a comparison, 
evolution of the MC1R gene in loaches does not appear to be 
affected by positive selection. 

As for different loach lineages, analyses based on branch 
models showed that both Cobitidae and Botiidae had much 


lower w values than the background lineages, with Cobitidae 
reaching a level of significant difference and Botiidae showing 
no significant difference. Since differences in w values reflect 
differences in the level of constraint (Rosenblum et al., 2004), 
the lower w values suggest that MC1R gene evolution in both 
Cobitidae and Botiidae might be strongly affected by purifying 
selection pressure. For Balitoridae and Nemacheilidae, 
however, both have higher w values than those of the 
background lineages, with the former at a significant level. This 
suggests that both families might have experienced a more 
relaxed selection pressure and relatively faster evolutionary rate 
than that of the background lineages. Among the four loach 
families, the botiid species exhibit the most regular color pattern, 
commonly with black vertical bars on the body sides. Species of 
Cobitidae have several types of color patterns, though typically 
show Gambetta pigment lines, as displayed in the genus 
Cobitis, which are defined as five color zones from the dorsal to 
lateral sides (Chen & Chen, 2005; Chen et al., 2015). For 
balitorid and nemacheilid species, the color patterns are 
extremely complicated and diverse, and include vertical bars, 
horizontal bands, spherical blotches, scattered spots, and 
irregular cloud-like patches. Consequently, we inferred that the 
relatively stable color pattern in Botiidae and Cobitidae might 
partially result from strong purifying selection pressure on the 
MC1R gene, whereas the diverse color patterns in Balitoridae 
and Nemacheilidae might be associated with relaxed selection 
pressure on this gene. 


Phylogenetic application of the MC1R gene 

As mentioned above, in the fish species analyzed so far, the 
MC1R gene is reported to be a conserved single copy and 
single exon protein-coding gene with a total length of more than 
800 bp, making it a good molecular marker for phylogenetic 
research (Li et al, 2007). Until this study, however, no 
phylogenetic application of this gene has been reported for fish. 
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In our research, partial MC1R gene sequences (591 bp) were 
adopted for establishing the phylogeny of Cobitoidea, which 
yielded a well-supported phylogenetic tree, except for Botiidae 
with low support value. The present topology was consistent 
with previous studies using multiple nuclear loci (Chen et al., 
2009; Liu et al., 2012; Mayden et al., 2009), and showed that 
Botiidae included two distantly related clades, corresponding to 
the two subfamilies Leptobotiinae and Botiinae proposed by 
Slechtová et al. (2006). The genetic distance between these 
two subfamilies was 0.092+0.013, higher than that between the 
two subfamilies of Balitoridae (0.073+0.010). Therefore, it is 
reasonable to divide Botiidae into these two subfamilies. Within 
Balitoridae, the three genera, Vanmanenia, Formosania and 
Pseudogastromyzon, were not monophyletic. In our previous 
study (Liu et al., 2012), the monophyly of Vanmanenia and 
Formosania could not be supported either. Further study is 
needed to confirm the phylogeny of these three genera using 
more species samples, molecular markers, and morphological 
characters. 

It is worth noting that the partial segment of the MC1R gene 
used in this study was easy to amplify using genomic DNA as a 
template with general PCR procedures, and easy to sequence. 
Therefore, this gene segment could be used as a good 
molecular marker for the phylogenetic study of fish species in 
the future. 
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